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Abstract 

We have generalized the MATHEMATICA function Apart from 1 to N dimension, the generalized 
function $Apart can decompose any linear dependent elements in V£ to irreducible ones. The 
elements in V* can be viewed as the corresponding propagators which involve loop momenta, and 
the decomposition will be useful when one tries to perform the loop calculations using the packages 
such as FlRE and REDUZE, which have implemented the integration by parts (IBP) identities 
and Lorentz invariance (LI) identities. A description on how to use this package, combined with 
Fire, FeynArts and FeynCalc packages, to do the one-loop calculations in double quarkonium 
production in e + e~ colliders is given, and the full source code for a specific process (e + e~ — > 
J ftp + T]c) is also available. 
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PROGRAM SUMMARY 



Title of program: $Apart 
Programming language: Mathematica 

Available from: http : //power . itp . ac . cn/~f engf eng/apart/ 

Computer: Any computer where the MATHEMATICA is running. 
Operating system: Any capable of running Mathematica. 
External routines: FeynCalc, FeynArts, FlRE 

Keywords: Next-to- Leading Order (NLO), Integrate By Parts (IBP), Apart 
Classification: 11.1 

Nature of physical problem: The traditional method to compute cross sections for a 
physical process in perturbative quantum field theory involves generating the ampli- 
tudes via Feynman diagrams and performing the dimensionally regularized loop inte- 
grals [![. Simplifications of the expressions are performed at the analytical level; there, 
an essential part is the reduction of these loop integrals to a small number of standard 
integrals. This step can be performed at the amplitude level for tensor integrals or, 
after contraction of Lorentz indices, at the level of interferences for scalar integrals. 
Considering the case of scalar integrals, integration by parts (IBP) identities 0, y] and 
Lorentz invariance (LI) identities [4j may be used for a systematic reduction to a set 
of independent integrals, called master integrals (MI). The standard reduction algo- 
rithm by Laporta [5|] defines an ordering for Feynman integrals, generates identities 
and solves the resulting system of linear equations. Alternative methods to exploit 
IBP and LI identities for reductions have been proposed see also fiol . and 

references therein. Public implementations of different reduction algorithms are avail- 
able with the comput er p rograms AIR [Hi, FIRE (Hi and Reduze Q. The usage 



of Fire[13J, Reduze|14|, etc. requires that the propagators must be decomposed to 
independent ones, for 1 dimension, there is a MATHEMATICA function Apart to do 
this, while for N dimension there is no such package yet, so we want to generalize the 
MATHEMATICA function Apart to $Apart in N dimension. 

Method of solution: We first prove that any linear dependent elements in V* can be 
decomposed into the summation of linear independent ones, the procedure of the proof 
gives us a method to perform the decomposition, $Apart is such an MATHEMATICA 
package that implements this method and generalizes the MATHEMATICA Apart func- 
tion from 1 to N dimension. 

Running time: Depends on the complexity of the system. 
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LONG WRITE-UP 



I. INTRODUCTIONS AND NOTATIONS 



A polynomial is a mathematical expression involving a sum of powers in one or more vari- 
ables multiplied by coefficients. A polynomial in one variable (i.e., a univariate polynomial) 
with constant coefficients is given by 



a n x n + • • • + a 2 x + a\x + ao ■ 



(1) 



The highest power in a univariate polynomial is called its order, or sometimes its degree, 
and for a polynomial with more than one variable, i.e. a multivariate polynomial, one needs 
to find the degree of each term by adding the exponents of each variable in the term, and the 
largest such degree is the degree of the multivariate polynomial. For a rational expression, 
i.e. an expression that is the ratio of two polynomials, we can work out its degree by taking 
the degree of the top (numerator) and subtracting the degree of the bottom (denominator), 
for example, 



Deg 



x 3 + Ax + 9 
x 5 + 2x + 1 



Deg[a; 3 + 4x + 9] - Deg[x 5 + 2x + l] = -2 



(2) 



It is well known that there is a function in M ATHEMATICA named Apart which rewrites a 
rational expression of a univariate polynomial as a sum of terms with minimal denominators 
and gives the partial fraction decomposition of the rational expression, for example, 



Apart 



(x — a)(x — b) 



+ 



1 



{b — a)x — b (a — b)(x — a) ' 



(3a) 



Apart 



(x — a)(x — b)(x — c) 



+ 



1 1 

(a — b)(b — c) (x — b) (a — c) {b — c) (x — c) 
1 



(a — b) (a — c)(x — a) 



(3b) 



There is no such similar function for bivariate polynomials yet; for example, Apart does 
not change the following form at all: 



Apart 



1 



x(x + a) (x + y + b) J x(x + a) (x + y + b) 



(4) 



We want to generalize this function to some specific rational expressions of n-variate poly- 
nomials; for example, we expect 



Apart 



1 



x(x + a)(x + y + b) 



{x,y} 



1 



1 



ax(x + y + b) a(x + a)(x + y + b) 



(5) 



First we introduce the notation. We consider the n-variate polynomial V* with degree 
less than or equal to 1. The n variables are denoted as {xj}i<j< n , and the linear space which 
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is spanned by n independent vectors {xi} over the coefficient field J 7 is denoted as V x . It is 
clear that 

V* = V x © T . (6) 

We call the k elements {e^ = Vi + fi}n<i<k) C V* are linear independent if and only if 
their projective parts in V x , i.e. {vi}, 1<i<k j C V^, are linear independent. 

Now we consider the following special terms generated from the rational operations on 
the polynomial V*: 

N 

JJe^, e^V^n^ZfyN^l , (7) 

i=l 

where Z is the integer set, and {^i} 1<i<N are generally not linear independent, i.e. reducible. 
We want to decompose them into a summation of linear independent, i.e. irreducible, ones: 

N Nj 

IT; >:/IK; ( § ) 

i=l j i=l 

where 1 < kji < N, 1 < Nj < N and for any fixed j, the Nj elements { e fc Ji } 1<i<Ar are linear 
independent, i.e. some elements in {ei} 1<i<N have been eliminated such that the remaining 
elements become linear independent. 

To give a proof of the decomposition in Eq. (jSJ), let us consider a special case: 

N 

F( ni ,--- t n N ) = l[e? , (9) 

i=l 

where any (N — 1) elements from {ei} 1<i<N are linear independent, but the N elements 
i e i}i<i<N are no ^' so there exists {fi}i<i<N with all ft ^ 0, such that 



X>; /• (io) 



t=l 



Note that the linear independency in V* is up to some constant / in the field J- '. 
We first look at Eq. with all = —1, i.e., 

F(-i,...,-i) = ii.-. — - (li) 

ei e 2 e^v-i e N 

• If / 7^ 0, we can write Eq. ( TTTT) as 

AT AT /TV \ AT TV TV 

n = nr M n r E 7 n ^' (12) 

i=i 1 ^ i=i * j \j=i j i=\ 1 j=i j i=i,i^j % 

since any (TV — 1) elements from {ei} 1<i<N are linear independent, the final expression 
is irreducible, and the desired decomposition. 
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• If / = 0, since all /, 7^ 0, without loss of generality, we take fi as an example: 

AT N / JV \ JV JV 7V 

n I = £i n I = I _l E/ , e ,) n I = _ E M JJ I. (13 ) 
f 4- e; ei - LJ - ej t\ \ U^ 1 / - LJ - /1 ef - L - L , e« 

i=l * 1 i=l 1 \ 3=2 / i=l j=2 J 1 1 i=2,ijtj 

We know any (iV — 2) elements from {ei} 2<i<iV combined with e\ are linear indepen- 
dent, so the final expression is also irreducible. 

So we get the decomposition for Eq. f JTTj) . Now considering the expression of Eq. ([9]) with 
all exponents of {e^} negative, 



JV N I 

G{n u --- ,n N ) = F(- ni ,--- ,-n N ) = Y[e^ = Y[ — , n t > , (14) 



i=l i=l 

n 1 

we can factorize out a term I [ — and perform the decomposition on it as follows: 

i=0 * 

• For the case when / 7^ 0, we have 

jv j jv ^ jv ^ jv ^ f N f- N 1 

i=l e * fc=i e fc i=l e « k =l e k \j=l J i=l,i^j e « 

JV JV 

3=1 j e j k=i,k& k 

where it is clear that we have decomposed the original term into N terms, and further- 
more that these terms have the same form as the original one except that one of the 
exponents rij decreases by 1 in each term, i.e., we get the following recursive relation: 

G (m, ■ ■ ■ , n N ) = 22 -J G ( n u ■■■ > n i ~ 1) • • ' , n N) (16) 
3=1 

and we can repeat the decomposition until one of rij decreases to 0. 

• When / = 0, 

jv 1 jv 1 n j jv j / jv ^ x jv x 

n ^ = n n r = n \ ~ n r 

i=l B< fc=i e fc i=l 1 k=i e k \ j=2 Jltl i=2,i^j 1 

N N 

j=2 - 71 1 e j k=l,i^j k 

and this is similar as the case / 7^ 0, the terms after decomposition have the same form 
as the original one, and one of the exponents rij decreases by 1 in each term except e\ 
whose exponent will increase by 1; the recursive relation is 

N fi 

G (m, ■ • • , n N ) = - y2 -yG (nt + 1, n 2 , ■ ■ ■ , rii - 1, • ■ ■ , n N ) (18) 
and we can repeat the decomposition until one of the rii{i > 2) decreases to 0. 
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So in each one of the two cases above, G (jii, • • • , n^) can be reduced to the summation of 
G (ni, n,2, • • • , rii = 0, • • • , njv) which can not be decomposed any more, i.e. it is irreducible 
and the desired result. 

If at least one exponent rii > in Eq. ([9]), without loss of generality, taking z = 1 and 
ni > 0, then the element e\ can be written as 



e i = ^ [f-Y.M > ( 19 ) 



we can substitute Eq. ([19]) into Eq. ([9]) to eliminate e\. 



N 
i=l 



n 



i=2 



W AT 



I e r = E/*IK fei ' ( 20 ) 



i=2 fc i=2 



and now the final expression only involves {ej}(2<i< n ) which are linear independent, and it 
is irreducible. 

To complete the proof, we will make the induction on N, i.e. the number of elements in 
{e{\. It is trivial that this is valid for N = 1, and now, assuming that it is also valid for 
N — 1, 2, 3, • • • , K, we want to prove that it is also valid for N — K + 1. 

If {ei} 1<i<K+1 are linear independent, i.e. irreducible, then there is no need for the de- 
composition; otherwise, there will be M + 1(M < K) elements from {e,} which are not 
linear independent, but any M elements are linear independent. Without loss of generality, 
we take these elements as {ej} 1<i<M+1 : 

K+l M+l K+l 

n ^ = n ^ n e ? ( 2i ) 

i=l i=l j=M+2 

Then according to the special case we have considered in Eq. ()9]), we have 

M+l Nj 

II < " Y. l U > ( 22 ) 

i=l j i=l 

where all Nj < M, so 

K+l I N 3 \ K+i / Nj K+l \ 

n = e a n n *=e* n ^ n * • < 23 > 

i=l V j 1=1 I m=M+2 j \i=l m=M+2 J 

Since Nj + (K - M) < K, i.e. the number of elements in [e kjl } 1<i<N . U { e m} M+2 < m <K+i in 
the right hand side(rhs) is less than N — K + 1, according to the assumptions, we have the 
following decomposition: 

Ni K+l K 

n^: n ^=E/in4'- < 24 > 

i=l m=M+2 fe t=l 
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with each term in the rhs irreducible, so we get the decomposition for N = K + 1: 



K+l 



(25) 



8=1 



k J 



i=l 



Since each term in r.h.s. is irreducible, the proof is done. 

The procedure also gives us a method to perform the decomposition. We will give an 
implementation in MATHEMATICA, i.e. the generalized Apart function: $Apart . 



II. AN IMPLEMENTATION IN MATHEMATICA 

The basic functions in the package are: 

• $Apart [expr , {x , y , z , . . . }] 

expr can be any form in Eq. ((7)), {x,y,z, . . .} are the corresponding n-variate poly- 
nomial variables, and $Apart [expr , {x , y , z , . . . }] will perform the decomposition on 
expr to give the irreducible forms, which are expressed with the function $ApartIR. 

• $ApartIR [expr , {x , y , z ,...}, {e t , e 2 e N > , {n.! , n 2 n N }] 
where expr is actually the product of e"*, i.e. 



expr 



N 

IK 



(26) 



1=1 



We preserve the {e 1 ,e 2 ,...,e N } and |rii ,n 2 n M } for later use, because they will be used 

as the input parameters for FlRE 13J. The irreducible form will be displayed as ||- ■ ■ ||, 
and $RemoveApart can be used to remove || in ||- ■ ■ ||. 

• $RemoveApart [expr] 

$RemoveApart is used to remove the HEAD in $Apart or $ApartIR, and is defined as 

$RemoveApart [expr_] : =expr/ . {$Apart [x_ , _] :>x, $ApartIR [x_ , ] : >x} . 

We can take Eq. (j3J) as the first concrete example: 



Apart 



1 



{x,y} 



1 




1 


x(b+x+y) 




(a+x)(b+x+y) 



x(x + a)(x + y + b) 

As another relatively complicated case, we take 

a 

expr 



(3a + b + c)(a + 26 + d) 3 (a + 46 + 9e) 



(27) 



(28) 
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If we only take a and 6 as the only polynomial variables, there are only two elements which 
can appear in the same irreducible expression, and we get the output as 



$Apart [expr , {a , b}] 



(2c - d) 



(3a+b+c)(a+2b+d) 3 



5(4c-9e) 



{3a+b+c)(a+2b+d) 2 



2c-lLi + 45e 



55(4c - 9e) 



1 



(3a+b+c)(a+2b+d) 



(2c- lld + 45e) 2 
121(4c-9e) 



(3a+b+c)(a+4b+9e) 



(2c- lld + 45e) ; 



(2c- lld + 45e) J 



2(2d- 9e) 



(a+2b+rf) 3 (a+4fe+9e) 



2(4c- 9e) 



i 



(a+2fe+d) 2 (a+46+9e) 



-2c + lid - 45e 



(2c- lld + 45e) 2 



22(4c - 9e) 



(a+2b+d)(a+4b+9e) 



(2c- lld + 45e) 3 
We can check the output with the original Eq. (1281) using the code 

dexpr = $Apart [expr , {a,b}] 
comp = expr - (dexpr//$RemoveApart) //Simplify 



(29) 



(30) 



The fact that comp gives zero indicates that the output dexpr is indeed identical with the 
original expr. 

If take c as a variable as well, we have 



$Apart [expr ,{a,b, c}] 



(3a + 6 + c)(a + 26 + d) 3 
+ (9e - 2d) 



+2 



(3a + b + c)(a + 26 + d) 3 (a + 46 + 9e) 
1 



(31) 



(3a + 6 + c)(a + 26 + d) 2 (a + 46 + 9e) 
More complicated examples can be found in Example/Examples .nb in the source code. 



"! 
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FIG. 1: A pentagon Feynman diagram for e + e — > 7* J/ift + rj c 

Now, let us look at how to apply $Apart to a specific feynman diagram from the process: 
e + e~ — > 7* — > J/ip + r] c which is shown in the Fig. [TJ After projecting the spin singlet 
and triplet for the charmonium r] c and J/ift respectively with the spin projectors 15|] and 
performing the DiracTrace on the fermion chains, we get the amplitude for this diagram as 



A 



mCACFegjei^Pi ((s - 4)sk 2 ml + k ■ p§ + k ■ p\ - (s - 2)k ■ p 3 k ■ p 4 ) 
3(D - 2)m c (s -A)sk 2 (k 2 + k ■ p 3 ) (4m2 + k 2 -2k ■ p 4 ) (2sm 2 + k 2 - k ■ p 3 - 2k ■ p 4 ) (k 2 - k ■ p 4 ) 



(32) 



s 



where p^ and p± are the momenta of J/i/j and r] c respectively, 7 and ip represent the polar- 
izations of 7* and J/ip respectively, k is the loop momentum, m c is the mass of the charm 
quark, and s is defined by 

s _ = (p 3 + p 4 ) 2 (33) 

Am? c Am? c 

If we take k 2 , k-p 3 and k -p^ as the 3-variate polynomial variables, then the amplitude A 
has the same form as Eq. ([7]) after expanding the numerator, so we can perform the $Apart 
operation on it: 

$ApartL4., {k 2 ,k -p 3 ,k • p 4 }] (34) 

The output looks like: 



4iC A C F eg* E ""^P3P4 
3(0 - 2)ml(s - 4)(s - 2)s 



fc 2 (fc 2 + k ■ p 3 ) (4m 2 + k 2 - 2k ■ p 4 ) 



m + 2s 



k 2 (k 2 + k ■ p 3 ) (2sm§ + k 2 - k ■ p 3 - 2k ■ p 4 ) || ° 



k 2 (4m 2 + k 2 -2k ■ p 4 ) (2sm 2 + k 2 - k ■ p 3 - 2k ■ p 4 ) 
1 



m + 2s 



+2(s - 2)s 
+2(s - 2) 
+ (4 - 2s) 

+ (4 - 2s) 



(fc 2 + fc ■ p 3 ) (4m 2 + fc 2 - 2fc • p 4 ) (fc 2 - fc • p 4 ) 
1 



fc 2 (4m 2 + fc 2 - 2fc ■ p 4 ) 
1 



(fc 2 + fc ■ p 3 ) (fc 2 - fc ■ p 4 ) 
1 



+ 4(s - 2) 



+ (2 - s) 



+ 2(s - 2)s 
1 



(fc 2 + fc ■ p 3 ) (4m 2 + fc 2 - 2fc • p 4 ) (2sm 2 + fc 2 - fc • p 3 - 2k ■ p 4 ) 
1 



fc 2 (2sm 2 + fc 2 - fc ■ p 3 - 2fc ■ pi) (fc 2 - fc ■ pi) 



(fc 2 + fc ■ p 3 ) (2sm| + fc 2 - fc • p 3 - 2fc ■ pi) 



+ (2 - s) 



fc 2 (fc 2 - fc • p 4 ) 



(4m| + fc 2 - 2fc ■ p 4 ) (fc 2 - fc ■ p 4 ) 



(2sm§ + fc 2 - fc ■ p 3 - 2fc ■ p 4 ) (fc 2 - fc ■ pi) 



(3r,) 



We can see that the 5-point integrals have been reduced to the 3-point or 2-point integrals, 
so the scalar integrals have been greatly simplified. Details of this example can be found in 
Process/FC-43 .nb in the source code, where 43 is the sequence number of the corresponding 
diagram which has been generated by FeynArts[16[ Il7 |. 



III. APPLICATION TO PHYSICAL LOOP CALCULATIONS 

The traditional method to compute cross sections for a physical process in perturba- 
tive quantum field theory involves generating the amplitudes via Feynman diagrams and 
performing the dimensionally regularized loop integrals Simplifications of the expres- 
sions are performed at the analytical level; there, an essential part is the reduction of these 
loop integrals to a small number of standard integrals. This step can be performed at the 
amplitude level for tensor integrals or, after contraction of Lorentz indices, at the level of 
interferences for scalar integrals. Considering the case of scalar integrals, integration by 
parts (IBP) identities @, 0] and Lorentz invariance (LI) identities [4fl may be used for a 
systematic reduction to a set of independent integrals, called master integrals (MI). The 
standard reduction algorithm by Laporta [1| defines an ordering for Feynman integrals, gen- 
erates identities and solves the resulting system of linear equations. Alternative methods to 



exploit IBP and LI identities for reductions have been proposed [6H9J; see also |10l . Illy and 
references therein. Public implementations of different reduction algorithms are available 
with the computer programs AIR [l2j|, FIRE (l3| and REDUZE jl4j . 

As for the one-loop calculations, there are many automatic tools available to achieve the 
general one-loop amplitude, such as FeynCalc[18} and FORMCALcjllJ, which are based 



on the traditional Passarino-Veltman|20H23j reduction of Feynman graphs, which can be 
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generated automatically(FeynArtsfl6l. [l7| or QGRAf|24|). In order to produce numerical 
results, tensor coefficient functions are calculated using LoopTools[19|. See also Refs. (2HI 



26] and the references therein. 

In the last few years, several groups have been working on the problem of constructing 
efficient and automatized methods for the computation of one-loop corrections for multi- 
particle processes. Many different interesting techniques have been proposed: these contain 



numerical and semi-numerical methods[27H30] . as well as analytic approaches |3 lH34| that 
make use of unitarity cuts to build next-to-leading order amplitudes b y g luing on-shell tree 
amplitudes (35l. [36| . For a recent review of existing methods, see Refs. |371 . |38| . 

In this section, we want to use the $Apart and the Fire(13T| package combined with 
FeynArtsEM \Vq and FeynCalc|19| to perform the one-loop calculations; here, Fey- 
NARTs(l6l. Il7| and FeynCalc[19| are used to generate the Feynman diagrams, and to 
perform the Dirac Trace respectively, the rest, such as tensor or scalar integral reductions, 
etc. will be handled by the $Apart and Fire[13| packages. We will concentrate on the 
next-to-leading-order corrections in a s to double quarkonium production in e + e~ colliders, 
the basic procedure can be summarized as follows. 



1. Use the FeynArts|16L |17| package to generate all Feynman diagrams for the partonic 



process: e + e~ —> 7* — > cc + cc. 

2. Use FEYNCALc(l9l| to perform the DiracTrace and SU(N) color matrix trace. 

3. Make an expansion in the relative momentum of the quark and the anti-quark in 
corresponding quarkonium and project out S-, P-, D-,- ■ - Waves. 

4. Use the $Apart to decompose the linear dependent propagators to independent ones. 

5. Use the FlREjj3i] package to reduce the general loop integrals to master integrals (MI). 

6. Process the final results, e.g. to asymptotically expand the amplitudes or calculate the 
cross section. 

Taking the process e + e~ — > J ftp + r] c as an example, for which there are large discrepancy 
between the NRQCD leading-order predictions and experimental data, and the important 
key step to resolve the discrepancies is that a large K factor of about 1.96 has been found 



in the next-to-leading-order corrections in a s |39l |40 . 

The calculations at leading order can be found in the directory Process/Tree, and the 
next-to-leading-order calculations can be separated into several parts: 

• Process/Tree/FC-RN .nb is used to calculate corrections from the counter-terms where 
the multiple renormalization is used. 

• Process/FC.nb is used to calculate the general loop corrections. 

• Process/FC-Nf .nb is used to calculate the corrections from the light quarks which 
are proportional to (Nf — 1). 

• Process/Total .nb will process the results generated from the above to give the nu- 
merical predictions or plots. 
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To compare with the results which are already present in other references, let us list some 
results which can be found in Process/Total .nb. 

The asymptotically expanded amplitude at s ^> 1, with s defined at Eq. (133]) . is 



A = *A (0) + ^U (1) + O (a 2 s ) 

A™ _ 1 
A® ~ 72 



39 hr s - 9(3 + 10 In 2) In s + 300 In— + 3(195 - 53 In 2) In 2 - 2tt z 



92 



ITi 



+— (-261ns + 301n2 + 9) + 
24 



(36) 



This result agrees with Eq. (6.4) of Ref. |4l| in which only the real part of asymptotic 
expansion is given, while our result also includes the imaginary part. 

We use the same input parameters as Ref. |!o[ to give the numerical results which are 



shown in Fig. [2] and Table HJ these are consistent with Refs. |39l |40 




15 20 
Q (GeV) 



FIG. 2: Cross section for e + e 



J/ip + rj c as function of the center-of-mass energy. The renor- 

1.5 GeV, the vertical line 



malization scale \i is set at half of the center-of-mass energy and m, 
corresponds to Q = 10.6GeV, solid line to oxo and dashed to onlo- 



m c (GeV) 


A» 




^Lo(fb) 


o-NLo(fb) 


O'NLo/o'LO 


1.5 


m c 


0.369 


16.09 


27.51 


1.710 


1.5 


2m c 


0.259 


7.94 


15.68 


1.975 


1.5 


Q/2 


0.211 


5.27 


11.14 


2.113 


1.4 


m c 


0.386 


19.28 


34.92 


1.811 


1.4 


2m c 


0.267 


9.19 


18.84 


2.050 


1.4 


Q/2 


0.211 


5.76 


12.61 


2.189 



TABLE I: Cross sections with different charm quark mass m c and renormalization scale /i, the 
input parameters are the same as Ref. [40| and Q = 10.6GeV. 



The method can be also used for calculations involving P-waves |42| 
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IV. SUMMARY 



We have introduced a generalized MATHEMATICA Apart function, which will perform 
the decomposition on any linear dependent elements in V* to reduce them to the irreducible 
form. The elements in V* can be viewed as the corresponding propagators which involve 
loop momenta, and the decomposition will be useful when one tries to perform the loop 
calculations using the packages such as the FlRE and the REDUZE, which have implemented 
the integration by parts (IBP) identities and Lorentz invariance (LI) identities. A description 
of how to use this package, combined with FlRE, FeynArts and FeynCalc packages, to 
do the one-loop calculations in double quarkonium production in e + e~ colliders is given, 
and the full source code for a specific process: e + e~ — > J/ift + rj c , from generating Feynman 
diagrams to the asymptotic expansion of amplitudes and numerical results, is also available. 
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